#核心论文代码验证
using Plots, ComplexPhasePortrait, ApproxFun, SpecialFunctions
gr();
u is in complex plane.
#u is in complex plane,for lanta_0,it has branch cuts
M0=0.3
phaseplot(-3..3, -3..3, z -> sqrt((1-z*M0)^2-z^2))
The fourier transform of the hard-wall boundary condtion on the centre body is simply 满足. To satisfy this boundary condition, we choose to write the solution of Bessel eqution as follows:
where and are the Bessel and Neumann functions of order m, and is the Hankel function of the first type.
参考:关于bessel的导数:https://www.boost.org/doc/libs/1_57_0/libs/math/doc/html/math_toolkit/bessel/bessel_derivatives.html http://functions.wolfram.com/Bessel-TypeFunctions/BesselJ/20/ShowAll.html
#首先让我们大致看一下bessel函数的曲线形式
#与https://dlmf.nist.gov/10.3一致
𝐽𝑚=(m,x) -> besselj(m, x)
𝑌𝑚=(m,x) -> bessely(m, x)
x = range(0.01,stop=10,length=1000)
plot(x, 𝐽𝑚.(0,x);label="J_0",ylims = (-1,1))
plot!(x, 𝐽𝑚.(1,x);label="J_1")
plot!(x, 𝑌𝑚.(0,x);label="Y_0")
plot!(x, 𝑌𝑚.(1,x);label="Y_1")
m=1
𝐽′𝑚=(m,x) -> 1/2*(besselj(m-1, x)-besselj(m+1, x))
𝑌′𝑚=(m,x) -> 1/2*(bessely(m-1, x)-bessely(m+1, x))
x = range(0.01,stop=10,length=1000)
plot(x, 𝐽′𝑚.(m,x);label="J'_1")
plot!(x, 𝑌′𝑚.(m,x);label="Y'_1",ylims = (-1,1))
f = z -> hankelh1(3,z)
# set up plotting grid
z = range(0.01 ,stop=100, length=1000)
plot(z, real(f.(z)))
#验证汉克函数
𝐻𝑚=(m,z) -> besselj(m, z)+im*bessely(m, z)
@show 𝐻𝑚(2, 1+im)
@show hankelh1(2,1+im)
phaseplot(-3..3, -3..3, z -> hankelh1(2,z))
#验证汉克函数在r趋近与无穷大时,叙述为正和负的区别
@show hankelh1(1,-100im)
@show hankelh1(1,+100im)
必须要指出的是,hankel函数用近似的方法只能针对与实数轴,对于过大的正虚数,imag》20im,则将不成立! 原因是,hankelh1(1, 20im)-->0,而besselj(1, 20im)-->∞&bessely(1, 20im)-->∞,因此造成很大的数值误差. 举例说明:
#验证hankel函数的导数
𝐻1𝑚=(m,z) -> hankelh1(m, z)
H1𝑚=(m,z) -> besselj(m, z)+im*bessely(m, z)
𝐻′1𝑚1=(m,z) ->hankelh1(m-1, z)-m/z*hankelh1(m, z)
H′1𝑚1=(m,z) ->besselj(m-1, z)+im*bessely(m-1, z)-m/z*(besselj(m, z)+im*bessely(m, z))
#H′1𝑚2=(m,z) ->-(besselj(m+1, z)+im*bessely(m+1, z))+m/z*(besselj(m, z)+im*bessely(m, z))
#𝐻′1𝑚2=(m,z) ->-hankelh1(m+1, z)+m/z*hankelh1(m, z)
@show besselj(1,2+20*im)
@show bessely(1,2+20*im)
@show 𝐻1𝑚(1,2+20*im)
@show H1𝑚(1,2+20*im)
@show besselj(1,2-20*im)
@show bessely(1,2-20*im)
@show 𝐻1𝑚(1,2-20*im)
@show H1𝑚(1,2-20*im)
zr=range(-40,stop=40,length=1000)
zi=ones(1000,1)*im*20
z=zr+zi
r_j = z -> real(besselj(1,z))
i_j = z -> imag(besselj(1,z))
r_y = z -> real(bessely(1,z))
i_y = z -> imag(bessely(1,z))
r_h=z -> real(hankelh1(1,z))
i_h=z -> imag(hankelh1(1,z))
#plot(zr,r_j.(z);label="r-j")
#plot!(zr,i_j.(z);label="i-j")
#plot!(zr,r_y.(z);label="r-y")
#plot!(zr,i_y.(z);label="i-y")
plot(zr,r_h.(z);label="T-r-h")
plot!(zr,i_h.(z);label="T-i-h")
#plot!(zr,r_j.(z)+i_y.(z);label="F-r-h")
#plot!(zr,i_j.(z)+r_y.(z);label="F-i-h")
M0=0.3
m=2
w=1+im
λ0=(u)->sqrt((1-u*M0)^2-u^2)
𝐻1𝑚=(u) -> besselj(m,u)+im*bessely(m,u)
𝐻2𝑚=(u) -> hankelh1(m,u)
𝐻1𝑚_λ0w=(u) -> besselj(m, λ0(u)*w)+im*bessely(m, λ0(u)*w)
𝐻2𝑚_λ0w=(u) -> hankelh1(m, λ0(u)*w)
𝐻′1𝑚1_λ0w=(u) ->hankelh1(m-1, λ0(u)*w)-m/(λ0(u)*w)*hankelh1(m, λ0(u)*w)
H′1𝑚1_λ0w=(u) ->besselj(m-1, λ0(u)*w)+im*bessely(m-1, λ0(u)*w)-m/(λ0(u)*w)*(besselj(m, λ0(u)*w)+im*bessely(m, λ0(u)*w))
phaseplot(-20..20, -40..40, 𝐻′1𝑚1_λ0w)
phaseplot(-20..20, -40..40, H′1𝑚1_λ0w)
m=2
w=1+im
λ0=(u)->sqrt((1-u*M0)^2-u^2)
𝐻1𝑚_λ0w=(u) -> besselj(m, λ0(u)*w)+im*bessely(m, λ0(u)*w)
𝐻2𝑚_λ0w=(u) -> hankelh1(m, λ0(u)*w)
phaseplot(-20..20, -20..20, 𝐻2𝑚_λ0w)
phaseplot!(-20..20, -20..20, 𝐻1𝑚_λ0w)
The fourier transform of the hard-wall boundary condtion on the centre body is simply 满足. To satisfy this boundary condition, we choose to write the solution of Bessel eqution as follows:
where and are the Bessel and Neumann functions of order m, and is the Hankel function of the first type.
那么问题转化为,在复平面范围内,避开branch cut的情况下,如何保证虚数部分始终为正数? 其中,,即均大于0
θ=range(0.0,stop=π/2-0.01,length=100)
plot(θ,tan.(θ))
wr = 0:1:100
wi = 0:1:100
f = z -> z
r = z -> abs(f(z))
θ = z -> angle(f(z))
plot( contourf(wr, wi, r.(wr' .+ im.*wi); title="abs"),
contourf(wr, wi, θ.(wr' .+ im.*wi); title="phase"))
M0=0 #M0=0,u±=±1
wr = 4
wi = 0
w=wr+ im*wi
phaseplot(-3..3, -3..3, z ->w*sqrt((1-z*M0)^2-z^2))
#u is in complex plane,for lanta_0,it has branch cuts
#λ0=sqrt((1-z*M0)^2-z^2)
#λ1=sqrt((C1-z*M0)^2-z^2)
M0=0 #M0=0,u±=±1
M1=0.5
C1=0.8
wr = 4
wi = 0
w=wr+ im*wi
xx = -10:0.01:10
yy = -10:0.01:10
f = z -> w*sqrt((1-z*M0)^2-z^2)
r = z -> abs(f(z))
θ = z -> angle(f(z))
image=z -> imag(f(z))
plot( contourf(xx, yy, r.(xx' .+ im.*yy); title="abs"),
contourf(xx, yy, θ.(xx' .+ im.*yy); title="phase"),
contourf(xx, yy, image.(xx' .+ im.*yy); title="image"))
#u is in complex plane,for lanta_0,it has branch cuts
#λ0=sqrt((1-z*M0)^2-z^2)
#λ1=sqrt((C1-z*M0)^2-z^2)
M0=0.3 #M0=0,u±=±1
M1=0.5
C1=0.8
wr = 4
wi = 4
w=wr+ im*wi
xx = -10:0.01:10
yy = -10:0.01:10
f = z -> w*sqrt((1-z*M0)^2-z^2)
#r = z -> abs(f(z))
#θ = z -> angle(f(z))
image=z -> imag(f(z)) #求f函数的虚数部分
tanε=wi/wr
u1ip=range(10,stop=-10,length=10)
u1rp=u1ip./-tanε+ones(10,1)*1/(M0+1)
u1im=range(10,stop=-10,length=10)
u1rm=-u1im./tanε+ones(10,1)*1/(M0-1)
u2ip=range(10,stop=-10,length=10)
u2rp=u2ip./-tanε+ones(10,1)*C1/(C1*M0+1)
u2im=range(10,stop=-10,length=10)
u2rm=-u2im./tanε+ones(10,1)*1/(C1*M0-1)
plot(contourf(xx, yy, image.(xx' .+ im.*yy); title="image"))
plot!(u1rp,u1ip,aspect_ratio=:equal;label="R1+0",arrow=1)
plot!(u1rm,u1im,aspect_ratio=:equal;label="R1-0",arrow=1)
#plot!(u2rp,u2ip,aspect_ratio=:equal;label="R2+0")
#plot!(u2rm,u2im,aspect_ratio=:equal;label="R2-0")
#phaseplot(-3..3, -3..3, z ->w*sqrt((1-z*M0)^2-z^2))
Finally,we can derive the Winer-Hopf equation:
M0=0.25
M1=0.5
D1=C1=1
m=2
w=4+4im
h=0.1
λ0=(u)->sqrt((1-u*M0)^2-u^2)
λ1=(u)->sqrt((C1-u*M1)^2-u^2)
tanε=imag(w)/real(w)
u1ip=range(10,stop=-10,length=10)
u1rp=u1ip./-tanε+ones(10,1)*1/(M0+1)
u1im=range(10,stop=-10,length=10)
u1rm=-u1im./tanε+ones(10,1)*1/(M0-1)
u2ip=range(10,stop=-10,length=10)
u2rp=u2ip./-tanε+ones(10,1)*C1/(C1*M1+1)
u2im=range(10,stop=-10,length=10)
u2rm=-u2im./tanε+ones(10,1)*1/(C1*M1-1)
#####
𝐽𝑚_λ1w=(u) -> besselj(m, λ1(u)*w)
𝑌𝑚_λ1w=(u) -> bessely(m, λ1(u)*w)
#𝐽𝑚_λ0w=(u) -> besselj(m, λ0(u)*w)
#𝑌𝑚_λ0w=(u) -> bessely(m, λ0(u)*w)
𝐽′𝑚_λ1w=(u) -> 1/2*(besselj(m-1, λ1(u)*w)-besselj(m+1, λ1(u)*w))
𝑌′𝑚_λ1w=(u) -> 1/2*(bessely(m-1, λ1(u)*w)-bessely(m+1, λ1(u)*w))
𝐽′𝑚_λ1wh=(u) -> 1/2*(besselj(m-1, λ1(u)*w*h)-besselj(m+1, λ1(u)*w*h))
𝑌′𝑚_λ1wh=(u) -> 1/2*(bessely(m-1, λ1(u)*w*h)-bessely(m+1, λ1(u)*w*h))
𝐻1𝑚1_λ0w=(u) -> hankelh1(m, λ0(u)*w)
𝐻′1𝑚1_λ0w=(u) ->hankelh1(m-1, λ0(u)*w)-m/(λ0(u)*w)*hankelh1(m, λ0(u)*w)
#𝐻′1𝑚2_λ0w=(u) ->-hankelh1(m+1, λ0(u)*w)+m/(λ0(u)*w)*hankelh1(m, λ0(u)*w)
K=(u)->D1*(1-u*M1)^2*λ0(u)*(𝑌′𝑚_λ1wh(u)*𝐽𝑚_λ1w(u)-𝐽′𝑚_λ1wh(u)*𝑌𝑚_λ1w(u))/(𝑌′𝑚_λ1wh(u)*𝐽′𝑚_λ1w(u)-𝐽′𝑚_λ1wh(u)*𝑌′𝑚_λ1w(u))-(1-u*M0)^2*λ1(u)*𝐻1𝑚1_λ0w(u)/𝐻′1𝑚1_λ0w(u)
phaseplot(-15..15,-10..10, K)
plot!(u1rp,u1ip,aspect_ratio=:equal;label="R1+0",arrow=1,legend=:bottomright)
plot!(u1rm,u1im,aspect_ratio=:equal;label="R1-0",arrow=1)
plot!(u2rp,u2ip,aspect_ratio=:equal;label="R2+0",arrow=1)
plot!(u2rm,u2im,aspect_ratio=:equal;label="R2-0",arrow=1)